Experiment & contact info

PIs: Rosalie Sears (searsr@ohsu.edu), Laura Heiser (heiserl@ohsu.edu)

Sample preparation: Zinab Doha (dohaz@ohsu.edu)

Library prep from single cells: Xi Li & Colin Daniel (danielc@ohsu.edu)

Analysis: Nick Calistri (calistri@ohsu.edu)

Sequencing performed by OHSU MPSSR

Set up

library(Matrix)
library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.8
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1

## Warning: package 'tibble' was built under R version 4.1.2

## Warning: package 'tidyr' was built under R version 4.1.2

## Warning: package 'readr' was built under R version 4.1.2

## Warning: package 'dplyr' was built under R version 4.1.2

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::expand() masks Matrix::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## x tidyr::pack()   masks Matrix::pack()
## x tidyr::unpack() masks Matrix::unpack()

library(Seurat)

## Warning: package 'Seurat' was built under R version 4.1.2

## Registered S3 method overwritten by 'spatstat.geom':
##   method     from
##   print.boxx cli

## Attaching SeuratObject

library(rliger)

## Loading required package: cowplot

## Loading required package: patchwork

## 
## Attaching package: 'patchwork'

## The following object is masked from 'package:cowplot':
## 
##     align_plots

library(SeuratWrappers)
library(ggalluvial)
library(harmony)

## Loading required package: Rcpp

so_list <- readRDS('analysis_files/s1_seurat_list.rds')

# Read in readme with HTO and tumor ID
lib_meta <- read_csv('original_data/README.csv')

## Rows: 78 Columns: 12
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (11): exa_dir, mpssr_seq_run_name, mpssr_library_name, library_suffix, c...
## lgl  (1): feature_barcodes
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Add the run_library variable
lib_meta <- lib_meta %>%
  mutate(run_library = paste0(mpssr_seq_run_name,
                              '_',
                              mpssr_library_name,
                              str_replace_na(library_suffix, replacement = '')))

lib_meta$mouse_id <- str_split(lib_meta$condition_string, pattern = '_', simplify = TRUE)[,1]

lib_meta$phenotype <- str_replace(lib_meta$mouse_id, pattern = "[:digit:]+", "")

tumor_histology_dict <- lib_meta %>%
  dplyr::select(condition_string, histology) %>%
  distinct()

Merge data and visualize QC metrics

so_merge <- merge(x = so_list[[1]], y = unlist(so_list)[-1])

## Warning in CheckDuplicateCellNames(object.list = objects): Some cell names are
## duplicated across objects provided. Renaming to enforce unique cell names.

# Add metadata: percent.mt, phenotype (from tumor ID)
so_merge[['percent.mt']] <- PercentageFeatureSet(so_merge, pattern = '^mt-')

so_merge[['percent.ribo']] <- PercentageFeatureSet(so_merge, pattern = '^Rp[sl][[:digit:]]')

so_merge[['phenotype']] <- str_replace(str_split(so_merge@meta.data$tumor,
                                                 pattern = '_',
                                                 simplify = TRUE)[,1],
                                       '[:digit:]+',
                                       '')

so_merge[['mouse']] <- str_split(so_merge@meta.data$tumor,
                    pattern = '_',
                    simplify = TRUE)[,1]

so_merge[['histology']] <- plyr::mapvalues(x = so_merge@meta.data$tumor,
                                           from = tumor_histology_dict$condition_string,
                                           to = tumor_histology_dict$histology)

## The following `from` values were not present in `x`: S8_T2, S9_T1, S9_T2, R4_T1, R4_T2, R4_T3, R5_T1, R5_T2, R5_T3, S8_T1, PD1_T1, PD2_T1, PD3_T1, PD4_T1, PD5_T1, PD6_T1

# Filter anything not assigned SP or SR status
table(so_merge@meta.data$histology)

## 
##  Doublet Negative       SP       SR 
##    74023    36538     4003    13632

so_merge <- subset(so_merge, subset = histology %in% c('SR', 'SP'))

# QC plots
FeatureScatter(so_merge, feature1 = 'nFeature_RNA', feature2 = 'nCount_RNA', group.by = 'library_id')+
  geom_vline(xintercept = 250)

FeatureScatter(so_merge, feature1 = 'percent.mt', feature2 = 'percent.ribo', group.by = 'library_id')

qc_features <- c('nFeature_RNA', 'nCount_RNA', 'nCount_HTO', 'percent.mt', 'percent.ribo')

for(i in qc_features){
  p1 <- VlnPlot(so_merge,
                features = i,
                group.by = 'library_id',
                pt.size = 0)+
    theme(legend.position = 'none')
  print(p1)
}

# What is our breakdown between library ID and HTO demultiplexed tumor ID
freq_table_og <- so_merge@meta.data %>%
  dplyr::select(library_id, tumor) %>%
  table()

pheatmap::pheatmap(freq_table_og,
                   display_numbers = freq_table_og,
                   main = 'Cell counts before QC')

og_counts <- tibble(tumor = colnames(freq_table_og),
                      count = colSums(freq_table_og),
                      fraction = colSums(freq_table_og)/sum(freq_table_og))

knitr::kable(og_counts)

so_noqc_meta <- so_merge@meta.data

QC plots for supplement

# Individual feature plots
## 720x540 output scaled to 2"

ggplot(so_merge@meta.data, aes(x = tumor, y = log10(nCount_RNA), fill = tumor))+
  geom_boxplot()+
  theme_bw()+
  ggtitle('Number of UMIs per cell')

ggplot(so_merge@meta.data, aes(x = tumor, y = nFeature_RNA, fill = tumor))+
  geom_boxplot()+
  theme_bw()+
  geom_hline(yintercept = 250, linetype = 'dashed')+
  ggtitle('Number of unique features (genes) per cell')

ggplot(so_merge@meta.data, aes(x = tumor, y = percent.mt, fill = tumor))+
  geom_boxplot()+
  theme_bw()+
  geom_hline(yintercept = 25, linetype = 'dashed')+
  ggtitle('Percent mitochondrial reads')

# removal criteria
## 720x540 scaled to 2"

ggplot(so_merge@meta.data, aes(x = tumor, fill = DF.classifications))+
  geom_bar(position = 'dodge')+
  theme_bw()

ggplot(so_merge@meta.data, aes(x = tumor, fill = nFeature_RNA > 250))+
  geom_bar(position = 'dodge')+
  theme_bw()

ggplot(so_merge@meta.data, aes(x = tumor, fill = percent.mt < 25))+
  geom_bar(position = 'dodge')+
  theme_bw()

Trim data based on QC metrics

# Subset based on nFeatureRNA > 250 and percent.mt < 25

FeatureScatter(so_merge,
               feature1 = 'nFeature_RNA',
               feature2 = 'percent.mt',
               group.by = 'library_id')+
  geom_vline(xintercept = 250)+
  geom_hline(yintercept = 25)

ggplot(so_merge@meta.data, aes(x = nFeature_RNA, y = percent.mt, color = tumor))+
  geom_point()+
  facet_grid(DF.classifications~library_id)+
  theme_bw()

meta_filtered <- so_merge@meta.data %>%
  filter(nFeature_RNA > 250) %>%
  filter(percent.mt < 25) %>%
  filter(DF.classifications == 'Singlet')

so_trimmed <- subset(so_merge,
                     cells = row.names(so_merge@meta.data)[! row.names(so_merge@meta.data) %in% row.names(meta_filtered)])

so_merge <- subset(so_merge,
                   cells = row.names(meta_filtered))

# Look at what is left
freq_table <- so_merge@meta.data %>%
  dplyr::select(library_id, tumor) %>%
  table()

pheatmap::pheatmap(freq_table,
                   display_numbers = freq_table,
                   main = 'Cell counts after QC')

new_counts <- tibble(tumor = colnames(freq_table),
                      count = colSums(freq_table),
                      fraction = colSums(freq_table)/sum(freq_table))

qc_counts <- full_join(x = og_counts,
                       y = new_counts,
                       by = 'tumor',
                       suffix = c('.og', '.qc'))

qc_counts <- qc_counts %>%
  dplyr::select(-fraction.og, -fraction.qc) %>%
  mutate(count.removed = count.og - count.qc) %>%
  mutate(fraction.remain = count.qc/count.og)

knitr::kable(qc_counts)

## Alluvial plot of library/mouse/tumor
meta_alluvial <- so_merge@meta.data %>%
  group_by(library_id, phenotype, mouse, tumor, HTO_classification.global, histology) %>%
  summarize(n_cells = n())

## `summarise()` has grouped output by 'library_id', 'phenotype', 'mouse',
## 'tumor', 'HTO_classification.global'. You can override using the `.groups`
## argument.

ggplot(meta_alluvial, aes(y = n_cells, axis1 = library_id, axis2 = mouse, axis3 = tumor, axis4 = histology, fill = histology))+
  geom_alluvium()+
  geom_stratum()+
  geom_label(stat = "stratum", aes(label = after_stat(stratum)))+
  theme_minimal()+
  scale_x_discrete(limits = c('library', 'mouse', 'tumor', 'histology'))

## Warning in to_lodes_form(data = data, axes = axis_ind, discern =
## params$discern): Some strata appear at multiple axes.

## Warning in to_lodes_form(data = data, axes = axis_ind, discern =
## params$discern): Some strata appear at multiple axes.

## Warning in to_lodes_form(data = data, axes = axis_ind, discern =
## params$discern): Some strata appear at multiple axes.

# Sample breakdown

sample_stats <- so_merge@meta.data %>%
  group_by(histology) %>%
  dplyr::select(histology, tumor) %>%
  table()

so_trimmed <- FindVariableFeatures(so_trimmed)
so_trimmed <- NormalizeData(so_trimmed)
so_trimmed <- ScaleData(so_trimmed)

## Centering and scaling data matrix

so_trimmed <- RunPCA(so_trimmed)

## PC_ 1 
## Positive:  Ppp1r14b, Hsp90aa1, Prdx1, Anxa2, S100a6, Anxa1, Hspd1, Mgst1, Dstn, Dbi 
##     Tubb4b, H2afz, Zfand5, Morf4l2, Klf4, Nhp2, Prdx5, Prdx6, Bst2, Mt1 
##     Ier3, S100a16, Ctsl, Cyc1, Rbp1, Ube2s, Mrpl34, Epcam, Tpm1, Tnfrsf12a 
## Negative:  Igkc, Cd74, Hba-a1, Hbb-bs, H2-Aa, Apoe, Igha, H2-Ab1, C1qa, Lyz2 
##     C1qb, Ctss, Tyrobp, Wfdc17, S100a9, C1qc, Ms4a7, Ccl12, S100a8, Retnlg 
##     Retnla, Wfdc21, Aif1, Stfa2l1, Lst1, Ccl6, Cstdc4, Ifitm6, Clec4d, Il1b 
## PC_ 2 
## Positive:  Col3a1, Col1a2, Adamts5, Fn1, Sparc, Abi3bp, Nid1, Rarres2, Col1a1, C4b 
##     Mfap5, Clec3b, Ogn, C3, Col14a1, Cd248, Prss23, Cxcl12, Tnxb, Igf1 
##     Svep1, Ctsk, Eln, Fndc1, Thbs2, Efemp1, Serpina3n, Igfbp7, Ebf1, Cd34 
## Negative:  Epcam, Krt8, Cldn3, Krt18, Stmn1, Cldn7, Cldn4, Cdh1, Pclaf, Car2 
##     Hist1h1b, Plet1, Hist1h3c, Birc5, Ltf, Tacstd2, Mki67, Rrm2, Rab25, Top2a 
##     Hmgb2, Krt19, S100a14, Fxyd3, Bspry, Kif11, Wfdc18, Elf5, Spc24, Cdca8 
## PC_ 3 
## Positive:  Nkg7, Xcl1, Cd7, Klra5, Hist1h1b, Il2rb, Ctsw, Hist1h3c, Ptprcap, Gzmb 
##     Top2a, Skap1, Pclaf, Cd226, Gimap4, Hist1h3e, Hist1h4d, Hist1h1a, Gm19585, Smc2 
##     Mki67, Fbxo5, Kif11, Kif15, Hist1h2bm, Sla2, Hist1h2af, Inpp4b, Hist2h2bb, Coro1a 
## Negative:  S100a14, Cldn3, Krt8, Cldn4, Krt18, Tacstd2, 2200002D01Rik, Clu, Epcam, Cldn7 
##     Krt19, Ltf, Fxyd3, Hspb1, Lgals7, Sprr2a3, Rab25, Wfdc18, Plet1, Krt7 
##     Asprv1, Klf5, Cystm1, Arg1, Gjb2, Gm41386, Stra6, Ovol1, Elf3, Klk10 
## PC_ 4 
## Positive:  Pcolce2, Tnxb, Cd34, Opcml, Clec3b, Efemp1, Dpp4, Pi16, Igfbp6, Ly6c1 
##     Efhd1, Itih5, Heg1, Cmah, Pla1a, Lrrn4cl, Cadm3, Ebf2, Ccl11, Sema3c 
##     Pamr1, Scara5, Ifi205, Slfn5, Tppp3, Dpt, Anxa3, Mfap5, Ace, Pcsk6 
## Negative:  Slc7a2, Ecrg4, Mme, Mgp, Edil3, Il17b, Gpm6b, Ptn, Matn2, Smoc1 
##     Ndufa4l2, Col11a1, Serpine2, Spp1, Mmp3, Ank, Enpp1, Cox4i2, Col8a1, Tnfrsf21 
##     Spon1, Nrxn1, Vcam1, Ptgfr, Gas1, Col4a2, Cx3cl1, Tgfb2, Thsd7a, Igfbp2 
## PC_ 5 
## Positive:  Rrm2, Birc5, Prc1, Top2a, Pclaf, Hist1h3c, Hist1h1b, Nusap1, Cdk1, Cdca8 
##     Cenpw, Smc2, Esco2, Sgo1, Cdca3, Kif11, Ube2c, Ccnb2, Ccna2, Klra5 
##     Stmn1, Gzmb, Pcolce2, Itih5, Hmmr, Racgap1, Tpx2, Xcl1, Spc24, Diaph3 
## Negative:  Cd14, Il1b, Cxcl2, Ccrl2, Plek, Nlrp3, Clec4e, Clec4d, S100a9, Lilrb4a 
##     Acod1, S100a8, Lilr4b, Samsn1, Ccr1, Il1rn, Ifitm1, G0s2, Wfdc17, Alox5ap 
##     Mxd1, Slc7a11, Slc15a3, Hdc, Fcer1g, Trem1, Stfa2l1, Tyrobp, Spi1, Hcar2

ElbowPlot(so_trimmed)

so_trimmed <- RunUMAP(so_trimmed, dims = 1:10)

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 22:08:42 UMAP embedding parameters a = 0.9922 b = 1.112

## 22:08:42 Read 3593 rows and found 10 numeric columns

## 22:08:42 Using Annoy for neighbor search, n_neighbors = 30

## 22:08:42 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 22:08:43 Writing NN index file to temp file C:\Users\Nick\AppData\Local\Temp\RtmpmcWjm4\filec805e945b6c
## 22:08:43 Searching Annoy index using 1 thread, search_k = 3000
## 22:08:44 Annoy recall = 100%
## 22:08:44 Commencing smooth kNN distance calibration using 1 thread
## 22:08:44 9 smooth knn distance failures
## 22:08:45 Initializing from normalized Laplacian + noise
## 22:08:45 Commencing optimization for 500 epochs, with 150984 positive edges
## 22:08:53 Optimization finished

DimPlot(so_trimmed, group.by = 'phenotype')+
  coord_equal()

FeaturePlot(so_trimmed,
            features = c('nCount_RNA', 'nFeature_RNA', 'percent.mt', 'Hba-a1'),
            coord.fixed = TRUE)

FeaturePlot(so_trimmed,
            features = c('percent.mt', 'Hba-a1', 'Ptprc', 'Epcam'),
            coord.fixed = TRUE)

so_trimmed <- NULL
so_list <- NULL
gc()

##            used  (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells  3555951 190.0    6325171   337.9    6325171   337.9
## Vcells 58623460 447.3 2153601279 16430.7 2677994074 20431.5

QC stats table for supplement

# Stats table

so_tumor_stats <- so_merge@meta.data %>%
  group_by(tumor) %>%
  summarize(nCells = n(),
            mean_nCount_RNA = mean(nCount_RNA),
            mean_nFeature_RNA = mean(nFeature_RNA),
            mean_percent.mt = mean(percent.mt)) 

so_tumor_stats %>%
  DT::datatable()

# Statistics used in results text

## Number of tumors
length(unique(so_merge@meta.data$tumor))

## [1] 11

## Number of mice
length(unique(so_merge@meta.data$mouse))

## [1] 6

## Number of unique libraries
length(unique(so_merge@meta.data$library_id))

## [1] 10

## Mean and range number of cells per tumor
range(so_tumor_stats$nCells)

## [1]  490 5995

sum(so_tumor_stats$nCells)

## [1] 14042

## Mean genes per cell
mean(so_merge@meta.data$nFeature_RNA)

## [1] 1198.225

mean(so_merge@meta.data$nCount_RNA)

## [1] 3394.431

Integration with rliger and harmony (runs both)

if(file.exists('analysis_files/s2_vehicle_integrated.rds')){
  print('Loading existing s2_vehicle_integrated.rds file')
  so_merge <- readRDS('analysis_files/s2_vehicle_integrated.rds')
}else{
  print('so_merge_rliger.rds file not found.')
  print('Processing iNMF integration, 200k cells ~= 1 hour')
  
  #iNMF integration with rliger
  so_merge <- NormalizeData(so_merge)
  so_merge <- FindVariableFeatures(so_merge)
  so_merge <- ScaleData(so_merge,
                        split.by = 'library_id',
                        do.center = FALSE)
  
  so_merge <- RunOptimizeALS(so_merge,
                             k = 50,
                             lambda = 5,
                             split.by = 'library_id',
                             nrep = 5)
  
  so_merge <- RunQuantileNorm(so_merge,
                              split.by = 'library_id')
  
  # UMAP
  so_merge <- RunUMAP(so_merge,
                      dims = 1:ncol(so_merge[['iNMF']]),
                      reduction = 'iNMF')  
  
  # Harmony integration
  so_merge <- RunPCA(so_merge)
  
  so_merge <- RunHarmony(so_merge, 
                       group.by.vars = 'library_id',
                       plot_convergence = TRUE,
                       max.iter.harmony = 20)
  
  
  # Save output
  saveRDS(so_merge, file = 'analysis_files/s2_vehicle_integrated.rds')
}

## [1] "Loading existing s2_vehicle_integrated.rds file"

UMAP assessment of iNMF integration

DimPlot(so_merge,
        label = TRUE,
        group.by = 'library_id')+
  coord_equal()

DimPlot(so_merge,
        label = TRUE,
        group.by = 'library_id')+
  coord_equal()+
  facet_wrap(~library_id)

DimPlot(so_merge,
        label = TRUE,
        group.by = 'histology')+
  coord_equal()

DimPlot(so_merge,
        label = TRUE,
        group.by = 'histology')+
  coord_equal()+
  facet_wrap(~histology)

DimPlot(so_merge,
        label = TRUE,
        group.by = 'tumor')+
  coord_equal()

DimPlot(so_merge,
        label = TRUE,
        group.by = 'tumor')+
  coord_equal()+
  facet_wrap(~tumor)

sessionInfo()

sessionInfo()

## R version 4.1.1 (2021-08-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] harmony_0.1.0        Rcpp_1.0.7           ggalluvial_0.12.3   
##  [4] SeuratWrappers_0.3.0 rliger_1.0.0         patchwork_1.1.1     
##  [7] cowplot_1.1.1        SeuratObject_4.0.4   Seurat_4.1.0        
## [10] forcats_0.5.1        stringr_1.4.0        dplyr_1.0.8         
## [13] purrr_0.3.4          readr_2.1.2          tidyr_1.2.0         
## [16] tibble_3.1.6         ggplot2_3.3.5        tidyverse_1.3.1     
## [19] Matrix_1.3-4        
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1          backports_1.3.0       plyr_1.8.6           
##   [4] igraph_1.2.7          lazyeval_0.2.2        splines_4.1.1        
##   [7] crosstalk_1.2.0       listenv_0.8.0         scattermore_0.8      
##  [10] digest_0.6.27         foreach_1.5.2         htmltools_0.5.2      
##  [13] fansi_1.0.2           magrittr_2.0.1        tensor_1.5           
##  [16] cluster_2.1.2         doParallel_1.0.17     ROCR_1.0-11          
##  [19] remotes_2.4.2         tzdb_0.2.0            globals_0.14.0       
##  [22] modelr_0.1.8          matrixStats_0.61.0    vroom_1.5.7          
##  [25] R.utils_2.11.0        spatstat.sparse_2.1-0 rmdformats_1.0.3     
##  [28] colorspace_2.0-3      rvest_1.0.2           ggrepel_0.9.1        
##  [31] haven_2.4.3           xfun_0.29             crayon_1.5.0         
##  [34] jsonlite_1.7.2        spatstat.data_2.1-2   survival_3.2-11      
##  [37] zoo_1.8-9             iterators_1.0.14      glue_1.6.1           
##  [40] polyclip_1.10-0       gtable_0.3.0          leiden_0.3.9         
##  [43] future.apply_1.8.1    abind_1.4-5           scales_1.1.1         
##  [46] pheatmap_1.0.12       DBI_1.1.2             spatstat.random_2.1-0
##  [49] miniUI_0.1.1.1        riverplot_0.10        viridisLite_0.4.0    
##  [52] xtable_1.8-4          reticulate_1.24       spatstat.core_2.4-0  
##  [55] rsvd_1.0.5            mclust_5.4.9          bit_4.0.4            
##  [58] DT_0.20               htmlwidgets_1.5.4     httr_1.4.2           
##  [61] FNN_1.1.3             RColorBrewer_1.1-2    ellipsis_0.3.2       
##  [64] ica_1.0-2             farver_2.1.0          R.methodsS3_1.8.1    
##  [67] pkgconfig_2.0.3       sass_0.4.0            uwot_0.1.11          
##  [70] dbplyr_2.1.1          deldir_1.0-6          utf8_1.2.2           
##  [73] labeling_0.4.2        tidyselect_1.1.2      rlang_1.0.1          
##  [76] reshape2_1.4.4        later_1.3.0           munsell_0.5.0        
##  [79] cellranger_1.1.0      tools_4.1.1           cli_3.1.1            
##  [82] generics_0.1.2        broom_0.7.12          ggridges_0.5.3       
##  [85] evaluate_0.15         fastmap_1.1.0         yaml_2.3.5           
##  [88] goftest_1.2-3         bit64_4.0.5           knitr_1.37           
##  [91] fs_1.5.2              fitdistrplus_1.1-6    RANN_2.6.1           
##  [94] pbapply_1.5-0         future_1.24.0         nlme_3.1-152         
##  [97] mime_0.12             R.oo_1.24.0           xml2_1.3.3           
## [100] hdf5r_1.3.5           compiler_4.1.1        rstudioapi_0.13      
## [103] plotly_4.10.0         png_0.1-7             spatstat.utils_2.3-0 
## [106] reprex_2.0.1          bslib_0.3.1           stringi_1.7.6        
## [109] highr_0.9             RSpectra_0.16-0       lattice_0.20-44      
## [112] vctrs_0.3.8           pillar_1.7.0          lifecycle_1.0.1      
## [115] BiocManager_1.30.16   spatstat.geom_2.3-2   lmtest_0.9-39        
## [118] jquerylib_0.1.4       RcppAnnoy_0.0.19      data.table_1.14.2    
## [121] irlba_2.3.5           httpuv_1.6.5          R6_2.5.1             
## [124] bookdown_0.24         promises_1.2.0.1      KernSmooth_2.23-20   
## [127] gridExtra_2.3         parallelly_1.30.0     codetools_0.2-18     
## [130] MASS_7.3-54           assertthat_0.2.1      withr_2.5.0          
## [133] sctransform_0.3.3     mgcv_1.8-36           parallel_4.1.1       
## [136] hms_1.1.1             grid_4.1.1            rpart_4.1-15         
## [139] rmarkdown_2.11        Rtsne_0.15            shiny_1.7.1          
## [142] lubridate_1.8.0